Reassessing Benzene Risks Using Internal 
Doses and Monte-Carlo Uncertainty Analysis 


Human cancer risks from benzene have been estimated from epidemiologicai data, with supporting 
evidence from animal bioassay data. This article reexamines the animal-based risk assessments 
using physiologically based pharmacokinetic (PBPK) models of benzene metabolism in animals 
and humans. Internal doses {total benzene metabolites) from oral gavage experiments in mice are 
well predicted by the PBPK model. Both the data and the PBPK model outputs are also well 
described by a simple nonlinear (Michaelis-Menten) regression model, as previously used by Bailer 
and Hoel [Metabolite-based internal doses used in risk assessment of benzene. Environ Health 
Perspect 82:177-184 (1989)1. Refitting the multistage model family to internal doses changes 
the maximum-likelihood estimate IMLEI dose-response curve for mice from linear-quadratic 
to purely cubic, so that low-dose risk estimates are smaller than in previous risk assessments. In 
contrast to Bailer and Heel's findings using interspecies dose conversion, the use of internal dose 
estimates for humans from a PBPK model reduces estimated human risks at low doses. 
Sensitivity analyses suggest that the finding of a nonlinear MLE dose-response curve at low 
doses is robust to changes in internal dose definitions and more consistent with epidemiological 
data than earlier risk models. A Monte-Carlo uncertainty analysis based an maximum-entropy 
probabilities and Bayesian conditioning is used to develop an entire probability distribution for the 
true but unknown dose-response function. This allows the probability of a positive low-dose 
slope to be quantified: it is about 10%. An upper 95% confidence limit on the low-dose slope of 
excess risk is also obtained directly from the posterior distribution and is similar to previous q^” 
values. This approach suggests that the excess risk due to benzene exposure may be 
nonexistent (or even negative) at sufficiently low doses. Two types of biological information about 
benzene effects—pharmacokinetic and hematotoxic—are examined to test the plausibility of 
this finding, A framework for incorporating causally relevant biological information into benzene 
risk assessment is introduced, and it is shown that both pharmacokinetic and hematotoxic models 
appear to be consistent with the hypothesis that sufficiently low concentrations of inhaled 
benzene do not create an excess risk. — Environ Health Perspect 104(Suppl 6): 1413-1429 (1996) 
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Introduction 

For risk managers and public health 
decision makers, the most important ques¬ 
tions about human cancer risks from expo¬ 
sures to low concentrations of benzene are 
a) Is there any excess risk at low exposure 
concentrations (c.g., below 1 ppm)? b) If 
so, how large is the excess risk, and how 
does it vary with changes in exposure 
concentration and timing? 


After more than 20 years of vigorous 
investigation using increasingly sophisti¬ 
cated techniques of molecular biology, 
computer and statistical risk models, and 
epidemiological data analyses, tile answers 
to these questions remain unknown. This 
paper briefly reviews how inferences about 
low-dose human cancer risks are tradition¬ 
ally drawn from animal bioassay data. It 


then suggests that a different approach is 
needed for benzene because benzene biol¬ 
ogy is inconsistent with the assumptions 
embedded in the usual approach. Three 
methods that may lead to more realistic 
and informative estimates of benzene risks 
are introduced as follows: 

• Use of Monte-Carlo uncertainty analysis 
to quantify the probability distribution 
of the entire dose-response function for 
benzene-induced tumors in male 
B6C3F1 mice. In contrast to previous' 
Monte-Carlo uncertainty analyses for 
benzene that treat uncertain model 
parameter values as random variables 
(/), here the true response probabilities 
in different experimental dose groups 
are treated as random variables. 

• Use of two species-specific physiologi¬ 
cally based pharmacokinetic (PBPK) 
models, one for mice and the other for 
humans, to quantify internal dose- 
response functions and to extrapolate 
risks from mice to humans. This 
approach leads to conclusions opposite 
to those previously reached using a 
PBPK model for mice and a traditional 
allometric scaling (2) approach to 
inrerspecies extrapolation. 

• A modeling framework for incorporat¬ 
ing additional biological knowledge— 
about cytotoxicity, cell kinetics, and 
genotoxicity, as well as the information 
about pharmacokinetics and metabo¬ 
lism contained in the PBPK model— 
into benzene dose—response models. 
The quanta! bioassay data displayed in 

Table 1 motivate and illustrate rhe meth¬ 
ods developed. These data are for the most 
sensitive species, strain, and sex found 
among many animal bioassay studies 
reviewed by Crump and Allen (3). They 
have been referenced by rhe Occupational 
Safety and Health Administration (OSHA) 
(4) in its regulation of benzene. 

Applied risk assessment attempts to 

draw inferences from data sets such as that 
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Table 1 . Oral gavage data for H6C3F1 male mice. 


x=Dose inmg/kg/day 
benzene 

n{x) = Number of animals 
exposed 

rM = Fraction of mice with 
tumors (squamous 
cell carcinomas) 


0 

25 

50 

100 

50 

50 

50 

50 

0 

0.08 

0.40 

0.74 


From Crump and Allen (3), Appendix B, reproduced by 
the U.S. EPA |5). 
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described in Table 1 to answer the follow¬ 
ing types of questions about risk: a) What 
is the excess risk at 25 mg/kg/day? 
Although the observed risk is 0.08, this 
observation contains sampling variability 
and ignores the information in the other 
data points. Additional analysis is needed 
to obtain a more informative answer. 

b) What is the excess risk in male mice 
exposed to 1 mg/kg/day? Answering this 
question requires extrapolation beyond the 
range of the data, and therefore requires 
some sort of risk model or assumptions. 

c) What is the excess risk in humans 
exposed via inhalation to 1 ppm for 
45 years, starting at age 20? This type of 
question is of most interest to risk man¬ 
agers. However, it cannot be answered 
purely on the basis of statistical analysis of 
the data. Instead, modeling assumptions 
are required to permit extrapolation across 
species, route of administration, and dose 
regimen. For applied work, it is essential to 
understand how robust risk conclusions are 
to plausible changes in the model assump¬ 
tions. If slightly more realistic assumptions 
lead to large changes in estimated risk, 
which turns out to be the case for benzene, 
then the uncertainty about risks may be 
driven more by model uncertainty than by 
data uncertainty. However, model uncer¬ 
tainty is more difficult to quantify and 
manage. It has often been disregarded in 
regulatory risk assessments of benzene. 

Traditional risk assessment, discussed 
next, infers from the data in Table I that 
there is a strong, positive dose—response 
relation for benzene that extends even to 
very low doses. However, traditional risk 
assessment incorporates assumptions (e.g., 
that risk is determined by total adminis¬ 
tered dose, or chat effects of exposure on 
cell proliferation and cytotoxicity may be 
ignored) that may not be appropriate for 
benzene. Reanalysis of the data in Table 1 
using more flexible methods indicates that 
there is no evidence for an increased risk at 
low exposure levels. 

Traditional Statistical Risk 
Assessment of Benzene 

A traditional regulatory risk assessment 
proceeds as follows: 

Step 1. A dose metric and a statistical 
risk model are chosen to relate the admin¬ 
istered dose history to tumor probability. 

Step 2. The parameters of the risk 
model are estimated from the available 
experimental data, e.g., by choosing para¬ 
meter values that maximize the likelihood 
of the data. 


Step 3. Uncertainty about the parame¬ 
ters is quantified using statistical theory. 
Uncertainty about the dose metric and the 
statistical risk model formula are usually 
either ignored or treated by making risk 
estimates for several different choices of 
dose metrics and models. 

Step 4. Risks are interpolated among 
dose levels and extrapolated beyond the 
conditions of the experimental data based 
on the chosen dose metric and risk model, 
e.g., to new species and new time patterns 
of dose administration. Specifically, the 
dose metric establishes which administered 
dose histories are considered to create equal 
risks, while the risk model quantifies how 
large these risks are. 

These steps will now be illustrated for 
the data in Table 1. 

Step 1. Selection of a Dose Metric 

A dose metric converts any specified 
detailed administered dose history (a time 
sequence of hourly or daily benzene intakes 
over a period of weeks to years) into a single 
corresponding number called the “dose,” 
typically plotted on the horizontal axis of 
the dose-response curve. Any two histories 
that are assigned the same number (dose) 
by the dose metric are thereby assumed to 
create the same risk, even if the time pat¬ 
tern or route of administration is different, 
and even if different species and strains are 
involved. An example of a dose metric that 
has been applied to the data in Tabic I (5) 
is mg-Ufecimes of benzene per unit of sur¬ 
face area of the exposed animal. For any 
dose history (e.g., 60 days of exposure to 25 
mg/kg/day), the corresponding dose num¬ 
ber is obtained by calculating the total 
quantity of benzene (e.g., 1500 mg), multi¬ 
plying it by the fraction of life exposed 
(e.g., 60 days/1000 days for a species and 
strain that lives 1000 days on average), and 
dividing by surface area (approximated as 
the two-thirds power of the animal’s body 
weight). In the remainder of this arricle, x 
will denote the dose variable defined by the 
selected dose metric (e.g., mg-lifetimes per 
unit surface area). 

The most common choice of dose 
metric for risk interpolation and extrapola¬ 
tion within a single species, strain, and sex 
is the area-under-curve (AUC) dose met¬ 
ric, which assigns the same dose (and 
hence the same risk) to any two dose his¬ 
tories that have the same integrated total 
dose (AUC for the administered dose 
history). For example, risk models in 
which equal ppm-hr of benzene or equal 
mg/kg of benzene create equal predicted 


risks independent of che detailed time 
pattern of dose administration are using 
AUC dose metrics. 

Step 2. Parameter Estimation via 
Maximum likelihood 

Once a dose metric has been selected, 
experimental data such as those in Table 1 
can be plotted to obtain an empirical dose— 
response function. A statistical model of 
the dose—response relation is typically rep¬ 
resented by a multivariate function (letters 
in bold type indicate vectors). 

P=f{x,q) 

where 

p - lifetime probability of tumor 

dose (computed via the dose metric 
selected in Step 1) 
q = a vector of model parameters 
f- assumed model formula. 

For example, the data in Table 1 might be 
assumed to result from the familiar multi¬ 
stage risk model, which has the form 

fixq) = 

l-e.\'p[-{q K) + q i x+q2X 2 + ... + q k x^)\, 
where 

qu Hz> ?*)- 

By varying q, the model dose—response 
function f(y.,q) may be made to approxi¬ 
mate the empirical dose-response func¬ 
tions obtained by plotting the fraction of 
animals with tumors against the dose, x. 
The quality of the approximation is quan¬ 
tified by a specific criterion, such as the 
likelihood criterion, which quantifies the 
likelihood of the empirical curve, as pre¬ 
dicted from the model curve; or the mean 
squared difference between che empirical 
and model values. Then q may be selected 
to optimize the selected criterion. This 
provides an estimate of q based on the 
observed experimental data, the assumed 
model, and the selected dose metric. 

As an example, the maximum-likeli¬ 
hood estimate (MLE) of q for the multi¬ 
stage model ivas computed for the data in 
Table 1, assuming an AUC dose metric of 
mg/kg/day of administered benzene. The 
value of q that maximizes the likelihood 
function was found via a numerical opti¬ 
mization computer program assuming that 
only parameter vectors satisfying the 
constraint q > 0 are ro be considered, as is 
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traditional in regulatory risk assessment. 
The MLE for q is 

2* = (?o, q !, q 2 )***(Q, 0.00145, 0.00013) 

(MLE parameter estimate). 

These parameter values maximize the like¬ 
lihood function 

L (q) ~ Pr (Data I q) 

= n x /(x,^)" (jt)rW [i-/u^)i [| - rWJ ” w , 

where n(x)r(x) denotes the number of ani¬ 
mals that developed tumors out of the n{x) 
animals in dose group x. They correspond 
to the MLE dose-response model 

p *{x) = 1 - exp [-(0.00145.v+ 0.00013X 2 )] 
(MLE risk model for mice). 

This model implies that, at very low doses, 

p*(x) = 0.00145* 

(MLE low-dose excess risk). 

Step 3. Uncertainty Analysis: 
Calculating qi* 

The quality of the MLE is often judged by 
uncertainty analysis. Traditionally, uncer¬ 
tainty analysis is based on construction of 
confidence intervals or joint confidence 
regions for the model parameters. An 
approximate confidence region for q may 
be constructed from the mathematical sta¬ 
tistics result that, under certain Technical 
regularity conditions, e.g., that L(q) is 
thrice differentiable and that q* lies in the 
interior of its feasible region, the distribu¬ 
tion of minus twice the log-likelihood 
ratio, namely 

-2 log [L(q)/Uq*)l 

asymptotically approaches a y 2 distribution 
with k degrees of freedom ((f). This 

asymptotic result can be used to construct 

approximate upper confidence limits 
(UCLs) for the true but unknown value of 
q from the MLE estimate q’. In particular, 
it can be used to estimate a 95% UCL for 
each component of q in the multistage 
model, including the linear term, qj. This 
provides the theoretical basis for calculat¬ 
ing q* in programs such as GLOBAL and 
TOXRISK. The same technique can be 
used to estimate approximate confidence 
bands for parameters in other risk models. 


The q* value calculated for the low-dose 
slope of the multistage dose—response func¬ 
tion for benzene is 

q\* = 0.00773 expected excess tumors 
per mg/kg/day of benzene. 

Equivalently, the estimated 95% upper 
confidence limit (UCL) on excess risk at very 
low doses is given by the linear function 

/>(*) =0.00773* 

(95% UCL on low-dose excess risk). 

Step 4. Extrapolating Risk from 
Oral (ravage Data for Mice to 
Inhalation Hazards for Humans 

To extrapolate from the oral gavage data 
for mice (Table 1) to risks posed by inhala¬ 
tion of low levels of benzene in humans, 
the U.S. Environmental Protection Agency 
(U.S. EPA) (5) selected mg-lifetimes of 
benzene per unit surface area as the equiva¬ 
lent dose metric for extrapolating risks 
across routes and species. They also used 
ppm-hr as a dose metric for risk interpola¬ 
tion and extrapolation within a species. 
The MLE for risk in any new situation 
(i.e., species, route, and/or time pattern of 
dosing) based on a risk assessment model 
/(*, q) and MLE parameter estimate q* 
assessed in a different situation is given by 
the risk extrapolation formula 

p*(x')=f(x',q*) 

where x denotes the dose corresponding to 
the new situation. For example, to estimate 
from the data in Table 1 the risk to a 
human male from occupational inhalation 
exposure to benzene, one would first calcu¬ 
late the “equivalent” mg/kg/day, x ', 
received by the worker (based on inhala¬ 
tion rate, body weight, and so forth) and 
then apply the dose-response formula 
derived from the mouse data 

/>*(*') = 1 -exp[— (0.00145*'+ 0.00013*' 2 )]. 

Thus, the risk created in a mouse by expo¬ 
sure to 1 mg/kg/day of benzene for a frac¬ 
tion /of its lifetime is assumed to be equal 
to die risk created in a human by continu¬ 
ous inhalation of the following concentra¬ 
tion of benzene: 

(20 m 3 inhaled by h umans/day) 
(1/70 kg for a human)(3-19 mg of 
benzene/m 3 of air inhaled per ppm 


of benzene in the inhaled air) 
(1// relative duration of lifetime 
human exposure) (70 kg per human/ 
mouse weight in kg) 1/3 , 

Using this factor for interspecics and inter- 
route dose extrapolation yielded a 95% 
UCL unit risk estimate for humans of 

0.119 excess tumors per ppm-ltfe- 
time of exposure to benzene (extrap¬ 
olated 95% UCL). 

The corresponding MLE for human 
risk due to inhalation of benzene was 
O.Oldl expected excess tumors per ppm- 
lifetitne of exposure (5). 

Critique of Traditional 
Benzene Risk Assessment 
and Need for a New Approach 

The preceding benzene risk assessment has 
three limitations. 

First, the assumption of an AUC dose 
metric for benzene is not realistic. Many 
experiments have shown that different time 
patterns of benzene dose administration 
with the same AUC produce very different 
profiles of benzert& metabolites ( 2 ) and 
very different hematotoxic effects (7). For 
example, 150 ppm of benzene inhaled for 

2 hr causes more than 7 rimes greater pro¬ 
duction of prephenylmercapturic acid (an 
indicator of benzene metabolism) than 50 
ppm of benzene inhaled for 6 hr (2). 
Similarly, exposing male CD-I mice via 
inhalation to 10 ppm of benzene for 
6 hr/day, 5 days/week for 10 weeks (3000 
ppm-ht of exposure) has no detectable 
impact on marrow cellularity or on colony¬ 
forming unit, granulocyte—macrophage 
(CFU-GM) scent cells in bone marrow, 
but 100 ppm for 6 ht/day for 5 days (also 
3000 ppm-hr of exposure) significantly 
depresses marrow colony-forming unit, 

spleen (CFU-S) and CFU-GM cells (7). 
Even more strikingly, while male NMRI 
mice continuously exposed to 21 ppm of 
benzene via inhalation for about a week 
show very significantly depressed marrow 
cellularity (cells/tibia) and CFU-GM con¬ 
tent per tibia, mice exposed to up to 14 ppm 
for up to 8 weeks — a much larger AUC 
dose-—show no significant changes in bone 
marrow cellularity or CFU-GM content 
(5). Othct AUC dose violations and anom¬ 
alies, such as the fact that exposure for 

3 days per week may have a larger impact 
on erythropoiesis than exposure for 5 days 
per week (E), have been documented for 
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f 


benzene metabolism, cytotoxicity, and 
genotoxidry. In light of these findings, die 
AUC dose metric cannot be considered 
biologically realistic For benzene. 

The second limitation is that the qf 
methodology may not produce productively 
useful results for benzene. Reasons include 
the following: a) The regularity conditions 
needed to establish the asymptotic result 
may not hold. For example, the usual regu¬ 
larity conditions require thac the true para¬ 
meter vector ^should lie in the interior of 
its set of possible values (so that the con¬ 
straint q> 0 is not binding). Yet, the MLE 
for q 0 is 0. b) Asymptotic convergence pro¬ 
vides no guarantees for real data set. Indeed, 
Monte-Carlo simulation shows that conver¬ 
gence is poor for some realistic-size data sets, 
while the nonnegativity constraints on 
model coefficients may lead to multi-modal 
distributions of q" values and to unstable 
MLE estimates of risk {10). The asymptotic 
^-square distribution does not always pro¬ 
vide a useful approximation to the actual 
distribution of q around its MLE estimate, 
as evaluated by Monte-Carlo analysis, 
r) Model uncertainty is ignored. The multi¬ 
stage mode! does not accurately describe 
some real data sets (11), This possibility is 
ignored in analyses of uncertainty using the 
qf methodology. Yet, such model uncer¬ 
tainty (e.g., about whether the chosen para¬ 
metric family of risk models contains the 
true dose-response relation) may overwhelm 
uncertainty due to statistical sampling vari¬ 
ability in determining total uncertainty 
about model-based risk predictions. 

These limitations imply that the com¬ 
puted value of q]* may contain very little 
information about the true plausible upper 
limit on the low-dose slope. A different 
uncertainty analysis methodology is needed 
to produce more useful results. Advances 
in uncertainty analysis methods, includ¬ 
ing the widespread use of Monte-Carlo 
uncertainty analysis and other computer¬ 
intensive statistical methods (12) now offer 
practical approaches for improving on 
earlier uncertainty analyses for benzene. 

The third limitation on benzene risk 
assessment is that tire human risk predictions 
made by the traditional risk assessment 
model based on animal data do not appear 
to agree with epidemiological data. For 
example, the preceding risk model predicts 
that the risk to a human exposed to x ppm- 
years of benzene via inhalation should exceed 

l-exp(-O.Oldx) 

(extrapolated lower bound 
dose-response model). 


(This is a lower bound because higher 
order terms are neglected.) Thus, a worker 
exposed to an average of 25 ppm-years per 
year for 40 years should have a lifetime 
tumor probability due to benzene exposure 
of at least 

l-exp[-(0.0l6)x(-1000 ppm-years/ 

70 years per lifetime)] = 0.20. 

This prediction appears to be incompatible 
with observed tumor rates in highly exposed 
human populations. For example, Turkish 
shoe and handbag workers exposed to an 
average of over 1000 ppm-years of benzene 
have a lifetime excess tumor probability 
that appears to be less than 0.01 ( 13 ), in 
contrast to the lower bound prediction 
from the above model. The 95% UCL unit 
risk estimates would tend to overpredict 
true risks even more. 

In summary, benzene risk assessment can 
be improved first by using a model of ben¬ 
zene risks that better reflects relevant ben¬ 
zene biology, The AUC dose metric is no 
longer plausible in light of current knowl¬ 
edge of benzene biology. The second way 
benzene risk assessment can be improved is 
by including model uncertainty in the 
uncertainty analysis. The possibility of 
model error must be addressed to accurately 
quantify benzene risks and uncertainties 
about them. 

Both of these improvements arc 
discussed in the following sections. The 
intended result is a set of benzene risk 
models that better agree with available 
human data and that more accurately 
predict human risks at low levels of 
benzene exposure. 

Methods For Improving 
Benzene Risk Assessment 

file following main technical ideas and 
methods may be used to construct new 
benzene risk models. 

Causal Decomposition of the Cancer 
Risk Process. Instead of trying to quantify 
the relation between dose and response 
probability directly, it is useful to decom¬ 
pose the causal relation between exposure 
history and tumor probability into biologi¬ 
cally meaningful causal links called “micro¬ 
relations” to quantify these links, and then 
to estimate the full dose-response relation 
by composing its constituent microrela¬ 
tions. As an example, the dose—response 
relation for benzene can be described as the 
composition of two micorelations, one 
linking administered dose to internal dose 


of benzene metabolites, the other linking 
internal dose to rumor probability. 

Model-free Curve Fitting. Instead of 
selecting the multistage model or another 
specific statistical risk model as the only 
model to be considered, it is potentially 
more realistic to express the cumulative 
tumor hazard a:; an unknown function of 
dose. If it is sufficiently smooth, this 
unknown function may be expanded as a 
convergent series (i.e., expressed as a 
weighted sum of orthonormal basis func¬ 
tions), and the coefficients of the first few 
terms may be estimated from data. This 
provides a data-driven, model-free approxi¬ 
mation to the true but unknown function 
(14). A similar but less sophisticated 
approach is to expand the unknown func¬ 
tion around zero as a MacLaurin power 
series. This has the advantage for risk 
analysis purposes that the power series 
model may then be interpreted as a multi¬ 
stage model with unconstrained coeffi¬ 
cients, i.e., coefficients not constrained to 
be nonnegative. Leaving the coefficients 
unconstrained may be more realistic than 
constraining them to be nonnegative if the 
effects of exposure on cell proliferation and 
cytotoxicity are considered. 

Maximum Entropy Bayesian Monte- 
Carlo Uncertainty Analysis. Instead of 
the q\* methodology’s asymptotic approach 
to uncertainty analysis, the following 
approach may be taken: a) Let the prior 
probability distributions for the true but 
unknown response probability in each dose 
group be uniformly distributed between 0 
and 1 (corresponding to a minimum of a 
priori assumptions). Use Bayes' Rule to 
condition this prior distribution on the 
observed data, i.e., on the number of ani¬ 
mals exposed and the number responding 
in each dose group. The result is a poster¬ 
ior distribution for the true but unknown 
response probability in each dose group. 
b) Sample many times from these distribu¬ 
tions (Monte-Carlo sampling) and pass a. 
model-free dose—response curve through 
each ser of sampled tumor probability 
values. In our implementation of this idea, 
polynomial recession was used to estimate 
the coefficients of the first few terms in the 
power series expansion of the cumulative 
tumor hazard function, considered as an 
unknown function of dose. The result is a 
Monte-Carlo distribution of the entire 
dose—response curve. If there were no 
uncertainty about the true response proba¬ 
bilities, they would uniquely determine a 
single dose—response curve, namely, the 
polynomial of lowest order that passe.; 
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through them. Because the true response 
probabilities are not completely deter¬ 
mined by experimental data due to sam¬ 
pling variability, a probability distribution 
of possible dose-response curves results. 
This distribution directly determines the 
upper 95% limit on the low-dose slope of 
the dose-response function, as well as a 
wealth of other information useful for 
characterizing uncertainty about the 
dose-response function. 

Biologically Based Computer Simula¬ 
tion Models. The causal decomposition 
framework suggests that the age-specific 
hazard rate for cancer (AML) induction by 
benzene depends on the administered dose 
history through two sets of variables: che 
internal dose of benzene metabolites (the 
output of a PBPK modet) on the one hand, 
and the product of hematotoxic and geno- 
toxic effects (if any) on the other. A non¬ 
linear feedback control model Is introduced 
to quantify the potential effects of benzene 
exposure on the granulocyte—macrophage 
(GM) lineage, specifically among prolifer¬ 
ating cells (e.g., CFU-GM cells) that may 
be involved in acute myeloid leukemia 
(AML) induction. Although this model has 
not yet been validated for benzene, its pre¬ 
dictions agree with empirical data in 
important respects, providing potential 
insights into why AUC dose metrics are 
inappropriate and how risk may vary with 
rbe time pattern of benzene administra¬ 
tion. Mode! results suggest that no single 
definition or measure of dose can ade¬ 
quately summarize the contributions to 
benzene-induced hematotoxicity and AML 
risk from different detailed histories of 
benzene administration. In general, che 
dose-time-response relation for benzene 
may be too complicated for a dose metric 
to exist. Instead, dynamic simulation mod¬ 
els must be used to calculate the effects 
of different administered dose histories. 
Details of these methods are discussed next. 

Maximum-entropy Bayesian 
Monte-Carlo Uncertainty 
Analysis 

Suppose that prior ro looking at the data, 
we wish ro make only a “minimum-com¬ 
mitment” assumption about the true value 
of the response probability fix) associated 
with dose x. This may be done by assign¬ 
ing to each p(x) a uniform probability 
density so that, a priori, p(x) is considered 
equally likely to have any of its possible 
values (between 0 and 1) (IS). Although 
stronger prior assumptions can be intro¬ 
duced by placing a priori constraints on die 


maximum-entropy joint prior density of 
the p(x) values—for example, by requiring 
the p(x) values to increase with x —such 
constraints may express requirements that, 
against a researcher's expectations, do not 
hold in biological reality. In the absence of 
definitive biological knowledge, it seems 
safer to avoid introducing such a priori 
assumptions, instead letting che shape of the 
p(x) curve be learned directly from the data. 
This may be done by conditioning the prior 
probability density function on the experi¬ 
mental data using Bayes’ Rule. Suppose the 
experimental data consist of the set 

D— {x, n(x), r(x) for M values of X | 
(experimental data) 

where M is the number of dose groups, 
then the revised a posteriori probability 
density for p(x) based on the observed data 
point [«(*), r(x )J (and assuming the maxi¬ 
mum-entropy, i.e., uniform, prior density 
U[0,1] with no assumed constraints inter¬ 
relating the p(x) values for different values 
of x ) is a beta distribution with parameters 
n(x)r(x ) + 1 and [ 1 - r( x) ] n(x) + 1. This 
has the mean value 

E [/>(*)! D] = E[p(x)\n(x), r(x)\ 

= [r(x)n(x) + \]/[n(x) +2] 
[posterior mean of/>(x , )l D], 

The beta posterior distribution for p(x) thus 
has a mean value that is dose to its observed 
sample value of r(x) when n(x) is large; 
moreover, the variance of the beta distribu¬ 
tion is a decreasing function of n(x) (16). 

Once che posterior distribution of each 
p(x) value has been determined, Af-vectors 
of possible response probabilities (one for 
each of the M dose groups) may be ran¬ 
domly sampled from these distributions. 
Each such AT-vector determines a corre¬ 
sponding dose—response curve, via any of 
several curve-fitting techniques. For illus¬ 
tration, we will use the relatively simplistic 
techniques (14) of polynomial regression. 

The mathematical derails are as follows: 
The tumor probability for an animal 
exposed to dose x may be expressed with¬ 
out loss of generality as: 

p{x) = 1 -exp[-A(x)], 

where h(x) is the cumulative tumor hazard 
associated with dose x. This is a mathemat¬ 
ical identity following from the definition 
of cumulative hazard. Next, the unknown 


cumulative hazard function is approxi¬ 
mated by the first few terms of a power 
series expansion around zero as: 

h(x) = q a + q x x+ q 2 x 2 +...+ q M - i* ' l#_l 

where M is the number of dose groups. 
For simplicity, the following exposition 
assumes that M = 4, as in Table 1. (The 
formulas hold for all M> 2.) The vector of 
cumulative hazards, denoted by 

h=(f> ] ,h 2 , h 3 , h 4 ), 

is determined from the values of the 
response probabilities (sampled from their 
appropriate (3 distributions) via the inverse 
transformation 

hj=-\n[\ -pj\, 

where pj denotes the sampled value of 
the rumor probability for dose groups'. 
This operation may be summarized by the 
equation 

A-~lti(l -p), 

if it is understood that che transform is 
applied componenc-by-component to the 
elements of p, the vector of response prob¬ 
abilities, to obtain the corresponding com¬ 
ponents of h. (This is consistent with the 
notation in several popular mathematical 
computing software packages, although it 
does not agree with standard vector algebra 
notation.) The coefficient vector 

q , = (qa><7n?2»?3) 

can now be solved for from h and the dif¬ 
ferent doses used in the experiment. The 
system of linear equations 

hi = q a + q\x { + q 2 x x 2 + /? 3 x, 2 

hi = <?o + + ?2*2 2 + ?3*2 2 

h ~ q u + + q 2 x 3 2 + qjX } 2 

hi = q a + q\X4 + q 2 x^ + q-sx? 

can be expressed compactly in vector- 
matrix notation as 

h^qX 

whereXis the matrix with (l, xj, Xj 2 , x 2 )^ 
as its j th column. (The equation /»= qX 
holds also for M > 4 dose groups, with the 
quantities q , and X extended in the 
obvious way.) X is known and h can be 
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Figure 1. Data flow for the uncertainty analysis. Probability distributions for the true response probability in each 
dose group (denoted by Pj for group j) are derived by conditioning maximum entropy (uniform) priors on observed 
experimental data, leading to beta distributions. [The cumulative density functions (CDFs), representing these beta 
distributions are shown in Figure 2 ] The CDF for each response probability' is transformed to a corresponding CDF 
for a cumulative hazard rate via the identity h,-= ln(l -/’■). The vector of cumulative hazard rates is used to solve 
for the vector of coefficients in a polynomial regression model via the formula ij=X- 1 h. Finally, the coefficients in 
q are used to compute the risk for each of several doses, x, via the formula p(xl= 1 -exp|-A(x)l = 
1 - expH<f[)+ <7]X+ q z i? +. . + q^’ |], 


estimated from the observed tumor inci¬ 
dence rates. The corresponding value of q 
is determined as follows: 

q = bX~*. 

The matrix inverse Jf' 1 can be found 
numerically; alternatively, the special struc¬ 
ture of X can be used to develop a simple, 
closed-form expression for the i/th dement 
of its inverse. 

Repeating the sampling of response 
probability M- vectors many times gives a 
Monte-Carlo sampling distribution of 
dose-response curves. Each sampled set of 
tumor probabilities, p, undergoes the 
following sequence of transformations: 

p-t-ln(l-p) = A-y bX~ ] = q -* p(x) 

= 1 - exp[- {q 0 + q x x + q 2 J- + q 3 x‘]. 

Thus, each p determines a corresponding 
dose-response curve expressed as a function 
of x at the right end of this sequence. 
Therefore, the beta distributions of the 
components of p imply a corresponding 
probability distribution for dose-response 
curves. This distribution of curves reflects 
the posterior uncertainty, after conditioning 
on the experimental data, about the true but 
unknown p{x) values. Figure 1 shows the 
data flow diagram for the entire algorithm. 

The Monte-Carlo distribution of dose- 
response curves allows a much richer char¬ 
acterization of uncertainties than does 
analysis. Not only can the value of the low- 
dose slope be determined such that 95% of 
all of the sampled dose—response curves 
have smaller slopes—the analog to the q' 
value in traditional analysis—but other 
detailed questions about the probable shape 
of the dose-response function can also be 
answered. For example, the probability 
distribution of the dose level at which a 
given level of risk (such as 1 £—06) is first 
exceeded can be determined from the 
Monte-Carlo sample of dose-response 
curves. The conditional probability distrib¬ 
ution of the risk at one dose level given 
assumptions about the risk at other dose 
levels can also be determined. 

Methods for a Biologically 
Based Causal Decomposition 
of Dose-Response Models 

The causal relation between the admin¬ 
istered dose history and the probability 
that a tumor develops by any specified time 
is thought to be mediated by several 
biological phenomena, including benzene 


metabolism and pharmacokinetics, ben¬ 
zene-induced cytotoxicity and genotoxicity 
in the bone marrow, and perhaps effects on 
the immune system. These phenomena are 
typically easier to observe experimentally 
than tumor risk itself. To discuss how they 
may be used to improve benzene risk assess¬ 
ment, it is useful to introduce the notation 

(a\\ b) = rime history of quantity a 

determined by the time history 
of quantity b, 

This natation is intended to imply that a 
and b are time-varying quantities, with the 

history of a up to and including any 
moment t being completely determined by 
the history of b up to the same moment. 
Other constants or parameters may have to 

be specified for the history {£(x), 0<T<r} 
to uniquely determine the history of {a(x), 
OSTS/l for every value of t such that 
T, where 0 is the time of the 
exposed animal’s birth and T is the time of 
its death. If so, these additional quantities 
may be listed after b following the “causal 
conditional” sign II in the above notation. 


For example, the dose-rime-response rela¬ 
tion between benzene exposure and AML 
risk may be denoted by 

(p\\x) 

to indicate that the probability of tumor by 
any time t is determined by the history of 
benzene exposures up to time t. Another 
suggested notation for a closely related 
concept is 

{*}->(/>}, 

read as “the history of x determines the his¬ 
tory of p .” By contrast, {/> 11.v} might be 
read as “the history of p char is determined 
by the history of x.” The curly brackets 
enclose individual time-varying quantities, 
or histories; thus, for example, \x } is an 
abbreviation of the more explicit notation 

[*(*), 0<r<T}. 

Now, suppose that AML risk depends 
on administered benzene only because of 
the formation of one or more benzene 
metabolites. Letting |y] denote the history 
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(i.e., the time course) of the vector of 
metabolites in different physiological com¬ 
partments over time, the causal situation 
may be diagrammed as 

In other words, {*} determines [y t and (y) 
determines fyij. Thus, the administered 
dose—time—response relation (y II x) may be 
decomposed as 

(?*}—> W) ‘ 

where * is rhe composition operator for 
composing consecutive mappings. This 
decomposition, based on internal dose 

of metabolites, may be expressed in the 
equivalent form 

(y>llx) = ( J oll_y)*(j»llx). 

The component (ytllar) corresponds to the 
input-output relation of a physiologically 
based pharmacokinetics (PBPK) model, i.e., 
it maps administered dose histories {*} into 
resulting time courses 5_y} of benzene and its 
metabolites in different physiological com¬ 
partments. The microrelation (p\\y) repre¬ 
sents an internal dose-response function. 

According to recent biological theories 
of carcinogenesis, most cancer risks can be 
expressed in terms of the following 
microrelations: 

a PBPK model, represented by the fink 

a cytotoxicity and cell kinetics model, 

a genotoxicity and cytogenetics effect 
model, 

a model of cell initiation, ]/V, pi -»| /]; 
a hazard rate model, |7, p} —> i h\, 
a ume-to-tumor model, \b] ->]/’}- 

In these links, the symbols are interpreted 
for benzene-induced AML as follows: \x] - 
administered dose history; {j} = metabo¬ 
lite history, giving the time Courses of con¬ 
centrations of benzene and benzene 
metabolites in different physiological com¬ 
partments (e.g., blood, bone marrow, fat, 
muscle, richly perfused tissue, poorly per¬ 
fused tissue, and various organ groups); 
[N] = time course of hematopoietic cell 
population sizes; ip! - history of cell transi¬ 
tion rate parameters (including birth, death, 
differentiation, and carcinogenic trans¬ 
formation rates); |7} = time course of a 
distinguished "initiated” (premalignant) 
cell population for AML; \h\ - cumulative 


hazard function, giving the age-specific or 
time-specific cumulative hazard for AML 
induction; and {/>} = history of tumor 
probability as a function of time. Vector 
quantities are indicated in boldface. 
Notation of the form {a,b}->{c] means 
that the history (cl is determined by the 
conjunction of histories {it} and {b\. Note 
that, as it is defined here, the quantity p(t) 
is greater than the probability of a 
detectable neoplasm being initiated by 
time t. To be detected, two further things 
must occur after a malignant cell forms; 
a) It must escape immune system surveil¬ 
lance and control; and b) it must progress 
to a detectable size. 

A Two-stage Stochastic Model 
of Cancer Induction 

The complete causal model for cancer 
induction in a two-stage stochastic model 
of carcinogenesis is composed from the 
preceding microrelations according to the 
following diagram; 

Here, the initiated population, 7, is inter¬ 
preted as a stochastic process (a birth- 
death process with immigration from the 
“normal” stem cell compartment), and 
only its probability law is determined by 
the data {p,/V}. The links in this causal 
chain may be quantified with the help of 
the following biomathematical model 
equations 

({*!-* W) [i] 

This link is quantified via the mathematical 
identity 

p(t) = \ - exp[-f>(r)l. 

Here, p{t) is the probability that at least 
one malignant cell is initiated by time t, 
Interpretively, if a randomly occurring 
event (formation of a malignant cell) 
occurs at a constant average rate of a times 
per year, then the probability that it has 
not occurred after t years is exp(— Ot/). The 
probability that it occurs at least once by 
time t is therefore 1 — exp(— ctf). The above 
identity generalizes this to the case of time- 
vatying occurrence rates, a(f)- The cumu¬ 
lative hazard b{t] is the integrated value of 
tt(r) (its AUC) up through time t. If Ot is 
constant, then h (f) = Ctr. 

{p 2 ,7}-*(/>} M 


This link is quantified via the formula 

dh(t)/dt= a(t) = y. 2 (t)I(t), 

where p 2 (f) is the transformation rate from 
initiated to malignant cells (per initiated 
cell per unit time) at time t. This formula 
has the common-sense interpretation that 
the expected number of malignant cells 
being formed at time L namely 0t( t \, is the 
product of the number of cells at risk of 
malignant transformation, 7(f), and the 
rate at which transformations occur, ( 12 (f)- 

(|i,Ar}-»{p 2 ,7} 13] 

This link is quantified via a stochastic 
differential equation with 

d£[I(t)]/dt= 

Wf)- 40- p 2 (f)U(f)+pi(fW(f)- 

Here, b(t) and d(t) are bitth and death 
rates and |ty(f) and p 2 (t) are transformation 
rates (all per cell per unit time); all four are 
components of the parameter vector p(f). 
Similarly, AT](t) refers to a component of 
N(t) (e.g., early myeloid hematopoietic 
progenitor cells). 

(y) -»{ji, iV} [4] 

This link describes the cytotoxic, cytoge¬ 
netic, and genotoxic effects of benzene 
metabolites on different cell populations. 
At present, sufficient in vivo experimental 
data and biological knowledge to reliably 
quantify this link are not available. Instead, 
therefore, the link 

1*} -»ip, atj 

will be quantified using curve-fitting tech¬ 
niques and experimental data. 

[ 3 ] 

This link can be quantified for benzene 
using any of several published PBPK mod¬ 
els of the dynamics of benzene metabolism, 
distribution, and elimination (77,/S). 

If each link can be quantified, then the 
links can be composed to obtain the total 
dose—time-response relation 

(14-»{/>» = (M -’tod) * (14->{p. N}) 
*({p,/V}->{p2,71) 

-►(«)* ««-> {/»)>• 
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This may be written more compactly as; 
0>ll*) = (/!l«‘(Hlp 2 ,/) 

*(p 2 ,/Hp,A0*( P ,Mij)*(jiij f ). 

However, this is only one of many possible 
decompositions. It is not the most useful 
one to use with existing benzene data since 
the component (p, Wily) is highly uncer¬ 
tain. Two alternative decompositions that 
can better use available data are 

-(p\\ h)* (h\\yL 2 ,I)* (p 2 ,/ II p, N) 

*(p. AT*). 

Different decompositions of the dose—time- 
response relation suggest different statistical 
strategies for quantifying it. For example, 
(/>llx), the undecornposed relation, 
may be estimated by fitting curves or 
models to experimental data on admin¬ 
istered doses and resulting tumors. This 
is the traditional approach illustrated 
already for the data in Table 1 using 
the U.S. EPA data ( 5 ) risk assessment 
as an example. 

(pll y) * (yll x) involves quantifying both 
the PBPK component, (y 11 .>;), based on 
experimental data on administered doses 
and resulting metabolite levels and also 

the internal dose-time-response rela¬ 
tion, (pll y), based on experimental data 
on metabolite levels and tumor rates, 
(pllA)*(AII>i a ,/)*(p 2> / lip, AO 

*(p, Afll x ) involves quantifying the 
relation (p, Affix) based on experimen¬ 
tal data on administered doses and cell 
level responses (cell proliferation rates, 
population sizes, cytotoxic impacts, and 
cytogenetic and genotoxic transforma¬ 
tion fares). For the remaining compo¬ 
nents, both ( p\\h) and (h\\\i 2 ,l) are 
obtained from the previously described 
mathematical equations. (p 2 , /ilp, N) 
involves identifying the cell population 
corresponding to / and the transforma¬ 
tion coriesponding to p 2 . 

Each of these strategies for decomposing 
and quantifying risk makes use of different 
data and has its own advantages. Each will 
be used subsequently in developing results 
on benzene risks. 

Quantifying Microrelations 
by Simulating Biological 
Processes: PBPK and Cyto¬ 
toxicity Models for Benzene 

To quantify a relation such as (yllx), it 
is usually necessary to build a dynamic 


simulation model that can accept time 
series histories {x} as inputs and calculate 
resulting output histories such as fy|. Such 
models typically take the form of systems 
of nonlinear ordinary or partial differential 
equations. In the simplest cases, the output 
history can be calculated from the input 
history' by iterating the following discrete¬ 
time approximation for small values of the 
time step, 

y{i + s)~y(t)+g\x(i) s y{t)]s 

(Euler’s method of numerical integration) 

where the instantaneous rate of change of y 
at time t is given by the formula 

dy/dt=g[x{t),y{t)]. 

The transition rate function g, indicated in 
bold because it is a vector lunction, deter¬ 
mines the time evolurion of y from its own 
current value and the concurrent value of 
the input, x. In practice, g"is obtained by 
biomathematical modeling of the specific 
biological system involved. 

Two examples of links for which dyna¬ 
mic models have been developed for benzene 
are PBPK models and cytotoxicity models. 

PBPK Models for Benzene 

The microrelation (yllx) is quantified hy 
PBPK models of benzene metabolism and 
distribution. Several such models have 
been published and are available for use 
(/,/7,/d). For risk assessment purposes, we 
prefer the human and animal PBPK models 
ofTravis et al. {17). These models quantify 
total benzene metabolites and do not track 
circulating metabolites. On the positive 
side, they have been extensively compared 
to and validated with human and animal 
experimental data. The details of the PBPK 
models are discussed by Travis et al. (17). 

Henderson et al. (19) present the fol¬ 
lowing three empirical y(x) data points for 
mice exposed to benzene via oral gavage, as 
in Table 1: >(12) = 10:y(40) = 25; and 
y(150) = 50. In each of these y(x) data 
points, the x value is the administered dose 
measured in mg of benzene per kg of ani¬ 
mal body weight and the y value is the 
resulting internal dose, defined as the total 
amount of metabolites formed in mg 
divided by the mouse body weight in kg. 
The interna! doses y(x) corresponding to 
the administered dose levels of x = 0, 25, 
50, and 100 mg/kg/day used in the NTP 
mouse bioassay experiment in Table 1 have 
not been directly measured. However, 
linear interpolation among these three 


published data points gives interpolated 

values of/(25) = 18.25; /(50) 30; and 
y(100) = 46. Alternatively, we can follow 
Bailer and Hoel (6) in assuming that the 
[x, y(x)\ points fall approximately on a 
nonlinear regression curve having the 
Michaelis-Menten form _y(x) = Vxl(K + x) 
where m = (V, K) is a parameter vector to 
be estimated from data. Least-squares esti¬ 
mation using the preceding three data 
points gives K‘ = 80.75, V* = 76.4, corre¬ 
sponding to y’(25) = 18.1, y *(50) = 29.2, 
and y*(100) = 42. These values are slightly 
less than the linearly interpolated values of 
y*( 25) = 18.25; y *(50) = 30; andy*(100) = 
46, The least-squares regression procedure 
smoothes out sampling variability and mea¬ 
surement error to the extent that the under¬ 
lying Michaelis-Menten model is correct. 
Given that the regression and interpolation 
values are close and that the former are 
more conservative, the Michaelis-Menten 
regression estimates could be used in the 
rest of the analysis. 

An alternative, more biologically based 
modeling approach is to use a benzene 
PBPK model to simulate the pharmacoki¬ 
netic and metabolic processes that convert 
administered doses to internal doses. Cur¬ 
rent PBPK models accomplish this simula¬ 
tion by representing the flow of benzene 
from compartment j to compartment / at 
time r(for nonmetabolizing compartments) 
by the equation 

fij(t) = k lj [C j (t)-C i (t)/P i ]. 

The rate of change in the concentration of 
benzene in compartment i at time l is thus 

dC i (t)7dt=Y. j [f^)-f ] ,(t% 

the flow into i from/minus the flow out of 
i into j. (If j is not adjacent to i, the flow 
between them is zero.) If benzene is metab ¬ 
olized in compartment i, such as liver or 
bone marrow, then dCj(t)ldt also includes 
a term [-K„ M ,C,-(r)/ P,V\K m + C ; (t))P,\. 
where V max and K m are the Michaelis- 
Menten parameters for the aqueous com¬ 
ponent of the compartment. Given the 
values of the pharmacokinetic parameters 
kij and P, and the metabolic parameters 
and K m , this system of differential 
equations can easily be integrated numeri¬ 
cally to solve for the time series | Q(r)). In 
current PBPK models for benzene, the 
metabolic parameters V max and k m are 
highly uncertain; they have been estimated 
by Medinsky et al. ( 18) and Travis et ai. 
(17) by seeking values that would make 
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model outputs match available data as well 
as possible. 

The bloodrtissue and blood:air parti¬ 
tion coefficients Pj for benzene are avail¬ 
able in the literature. Medinsky et al. (IS) 
provide values of bloodttissue and blood:air 
partition coefficients P, for i = liver, fat, 
poorly perfused tissues, richly perfused tis¬ 
sues, and air. Travis et al. (17) identified 
ranges of P, values from the literature for 
i = fat, liver, muscle, organs (corresponding 
to richly perfused tissues) and air and then 
chose values within these ranges that 
caused the PBPK model predictions to best 
fit available data. The flow rate parame¬ 
ters kjj are determined by PBPK theory 
as follows: 

• If the flow is from j = arterial blood 
to i= a tissue compartment, then 
k;j = Qj / V„ the blood flow rate (in liters 
per minute) through compartment i 
divided by the volume (in liters) of 
compartment i. Values of the physio¬ 
logical constants Qj and V t have been 
estimated for rats, mice, humans, and 
Other species (20). 

• If the flow is from j = alveolar air to 
/ = arterial blood (or vice versa), then 
the value of kg is assumed to be so large 
that the flow/7(f) = k,j[Cj(t) — C£t)tN\ 
can be modeled as adjusting to equilib¬ 
rium instantaneously. (N denotes the 
bloodrair partition coefficient.) Instead 
of having to estimate kjj, therefore, 
the equilibrium condition C a i v (i) = 
C m (t)/N is assumed. 

• The total flow of benzene into mixed 
venous blood from tissues per unit time 
is the sum over all tissue groups, i, of 
Q,Cj(r)/jP„ while the total flow of blood 
into the venous side of the circulation 
per unit time is just Q, the circulation 
flux (1/min), which is equal to the total 
cardiac output rate. If it is assumed that 
the volume of the venous system is neg¬ 
ligible, then the concentration of ben¬ 
zene in the venous blood will be the 
sum over i of (Qj!Q)C,(t)/Pi, i,e„ it 
will be the average concentration of 
benzene in blood entering the mixed 
venous blood. This determines the 
(approximate) value of C„ rn (t), the con¬ 
centration of benzene in mixed venous 
blood at rime f. 

In addition to these flow rate and equi¬ 
libration assumptions, the following flow 
balance constraint is also assumed: 

C m (t)(Q alv /N+ Q)=[CLiA„t,U) 


implying that 

Qn(r) = [CLtAnhit )+ £, Q/qM/fy / 

(O b/N+Q). 

This equates the total inflow of benzene 
per unit time into the alveolar-air-and- 
arterial-blood compartment to the total 
outflow of benzene per unit time from that 
compartment. 

These assumptions and parameter 
estimates allow the time series of benzene 
concentrations in each compartment to be 
determined from the inhalation exposure 
time series {C^ft)!. If the dose adminis¬ 
tration route is oral gavage or injection (sc 
or ip), the PBPK model must be extended 
to allow for transfer of the administered 
bolus dose into the rest of the system. For 
oral gavage, the simplest expedient is to 
assume a first-order gastric absorprion con¬ 
stant (transfer rate) from the GI tract to 
the liver. Both Medinsky et al. (IS) and 
Travis et al. (17) make this assumption 
and both estimate similar rate constants in 
mice—0.032/min according to Travis 
et al. and 0.05/min according to Medinsky 
et al. [However, the gastrointestinal trans¬ 
fer rate constant estimated for rats is 0.0042 
according to Medinsky et al. (IS) and 0.01 
according to Travis et al. (17). 1 

We have reimplemented the Travis et 
al. PBPK model for benzene (17), incorpo¬ 
rating the preceding modeling assumptions 
and parameter values. This model has been 
extensively compared againsc both human 
and animal data on time courses of ben¬ 
zene concentrations in expired air, blood, 
and bone marrow during and Following 
exposures to benzene via inhalation, oral 
gavage, and injection. It provides a useful 
lit to nearly all the available data sets, as 
discussed in detail by Travis et al. (17). The 
resulting simulation model for metabolism 
of oral gavage doses in mice predicts values 
for internal doses slightly greater than the 
ones from the Mlchaelts-Menten regression 
model. Table 2 compares the predictions 
from the PBPK and Michaelis-Menten 
internal dose models. 

The Michaelis-Menten model in Table 2 
was fit to the three empirical data points 
from which the parameter values K* = 
80.75 and l 7 * = 76-4 were estimated. A bet¬ 
ter fit to the PBPK model output is given by 
the revised Michaelis-Menten model 
y*(x) - 77.5x1(74.33 + x). This curve pro¬ 
vides a compact approximation of the PBPK 
results. At the National Toxicology Program 
(NTP) dose levels of x = 25, 50, and 


75 mg/kg, it predicts internal doses of 
J*(25) = 19.5, y*(50) = 31.2, and y*(100) = 
45-5, compared to the PBPK model predic¬ 
tions of y*(25) = 19.0, ^*(50) = 30.1, and 
J,*(1Q0) = 44.6. Agreement at other poincs is 
even closer. Thus, y*(x) = 77.5x7(74.33 + x) 
provides a useful analytic approximation 
to the PBPK model for predicting total 
metabolites formed from oral gavage doses. 

Cytotoxicity and Cell Kinetics 

Models for Benzene 

The link (iVllx) may be quantified by 
biomathematica! simulation models of 
granulopoiesis and hematotoxicity. Such 
models have been previously developed 
and tested against experimental data in 
mice for benzene (21) and in dogs for the 
immunosuppressive and myelotoxic agent 
cyclophosphamide (CP) (22). Although 
such models are relatively immature and 
exploratory compared to PBPK models, 
they do synthesize large amounts of rele¬ 
vant biological experimental data and may 
prove useful in future efforcs to improve 
benzene dose—time-response modeling. 

A computer simulation model of 
benzene-induced hematotoxicity is cur¬ 
rently being developed. While it would be 
premature to rely on it as a basis for ben¬ 
zene risk assessment, it offers potential 
insights that may be useful enough to jus¬ 
tify a brief outline of the model’s main fea¬ 
tures, It is based primarily on the model of 
Steinbach et al. (22) updated with benzene 
cytotoxicity parameters from Scheding et 
al. (21). We have validated the model with 
published human clinical data for CP 
(23). The main model consists of a set of 
compartments representing cell popula¬ 
tions, linked by flow equations with flow 
rate parameters that change over time in 
response to changes in the sizes of the cell 
populations. Several nonlineat feedback 
control laws operate to restore the cell 
population sizes to their stable equilibrium 

Table 2. Mouse internal doses of benzene metabolites 
(mg/kg) produced from different oral gavage doses, as 
predicted by various models. 


Dose Observed Predicted y\x) values 


level x 

values 

Michaelis-Menten 

PBPK 

1 

- 

0.93 

1 

5 

- 

4.5 

4.9 

10 

- 

8.4 

g .2 

12 

10 

9.9 

10.7 

25 

- 

18.1 

19.0 

40 

25 

25.3 

262 

50 

- 

29.2 

30.1 

ICQ 

- 

42.3 

44.6 

150 

50 

49.7 

54.5 
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values following any small perturbation, 
e.g., due to Inhalation of benzene. These 
feedback control laws are not powerful 
enough to maintain or restore equilibrium 
cell population levels if sufficiently intense 
exposure continues, however. Even after 
cessation of dosing, the system may not 
return to equilibrium for many weeks. 

The compartments of the cell kinetics 
and hematotoxicity model are defined and 
abbreviated as follows: 

• S = early myeloid stem cells and 
hematopoietic progenitor cells (H PCs). 
These cells are further subdivided into 
two subpopulations: 

cyc!ing_HPCs = actively cycling 
HPCs (i.e., those participating in the 
mitotic cycle); 

resting_HPCs = resring or dormant 
HPCs (not actively cycling). 

• CFU = granulopoietic committed stem 
cells. Biologically, these correspond 
roughly to CFU-GM cells, i.e., to early 
myeloid cells that can form GM 
colonies. To model accurately both the 
mean and the variance of transit times 
through this compartment, it is neces¬ 
sary to account for the age structure of 
the CFU-GM population. The partial 
differential equation (PDE) describing 
the time evolution of the age-structured 
population is approximated by dividing 
the compartment into 10 fictitious sub- 
compartments representing cells of dif¬ 
ferent ages. The age/maturity structure 
of this population also matters in deter¬ 
mining the dynamic response of the 
hematopoiecic syscem to CP-induced 
stresses. Subdividing it into 10 subcom¬ 
partments allows shifts in its age com¬ 
position to be approximated. 

■ P = proliferative pool, consisting 
roughly of myeloblasts, promyelocytes, 
and myelocytes. This aggregate com¬ 
partment is also subdivided into 10 
subcompartments to simulate the 
changing age structure of the prolifera¬ 
tive subpopulation over rime. 

• M = maturation pool, in which cells 
finish maturing and lose their remaining 
proliferation and differentiation capaci¬ 
ties, thus becoming no longer suscepti¬ 
ble to carcinogenic transformations. 

• R = bone marrow reserve of mature 
granulocytes. Granulocytes are released 
from this compartment to peripheral 
blood as needed to replace losses due to 
normal cell senescence and death or to 
emergencies such as bleeding. 

• B - mature granulocytes in circulating 
blood. In comparing model predictions 


to experimental data, this compartment 
represents an approximation to white 
blood cells (WBCs) (although lym¬ 
phoid cells are not modeled). 

These cell populations ate the compo¬ 
nents of the vector N in the microrelation 
(JVIl ac). Other aspects of granulopoiesis, 
such as the shift of some hematopoiesis to 
the spleen under stress, deterioration of the 
stromal microenvironment of the bone 
marrow over time, or migration of stem 
cells from the marrow inco the peripheral 
blood and back are ignored in this model 
since their significance for the dosing 
scenarios considered is expected to be small. 

The flow equations linking the com¬ 
partments of the model (22) express three 
key ideas: a) the flux of living cells out of 
each compartment enters its downstream 
neighbor; b) the flux of dead cells from 
each proliferative compartment is propor¬ 
tional to the number of cells currently in 
the compartment and to the concentration 
of administered benzene weighted by a 
compartment-specific cytotoxicity parame¬ 
ter; and r) die magnitude of the flux out of 
a compartment is directly proportional to 
the compartment’s current size (i.e., to the 
number of cells in it) and inversely propor¬ 
tional to the transit time of cells through 
the compartment. 

The number of cells in any of the com¬ 
partments at any time is found from the 
flow equations by incrementing the imme¬ 
diately previous number to reflect inflows 
from the upstream compartment and 
amplifications dne to cell divisions (for the 
early, proliferative compartments); and by 
decrementing the result to reflect outflows 
to downstream compartments and losses 
due to cell death. Cell death rates due to 
benzene metabolites (per susceptible cell 
per unit time) are assumed to depend on 
administered concentration of benzene, as 
discussed by Schcding et al. (2.1). 

The values of key flow rate parameters 
are determined via nonlinear feedback con¬ 
trol equations from the sizes of various cell 
populations. The hematopoietic system 
exerts several forms of feedback control to 
regulate its own behavior and compensate 
for cytotoxic and other stresses (21). Based 
on several decades of experimental data in 
animal models, the following five specific 
quantities are controlled by feedback 
control loops in the model. 

HPC_recruitment_rate = fractional 
recruitment rate of resting HPCs, defined 
as the fraction of resting HPCs recruited 
into active cycling per unit time. This rate 
is controlled not only by the concurrent 


numbers of resring and cycling stem cells, 
but also by the sizes of the downstream cell 
populations in the CFU-GM and subse¬ 
quent compartments up through and 
including the bone marrow reserve. If the 
more mature cell populations depart from 
their steady-state levels in response to 
hematopoietic stresses due to dosing, the 
recruitment and production of stem cells 
counter-adjusts to help move the system 
back toward its equilibrium levels of 
these populations. 

HPC__differentiation_fraction = The 
fraction of HPCs thar differentiate (rather 
than self-renewing). A simplified model of 
the cell cycle Is used for stem cells, indicat¬ 
ing the fraction (A/S) that are actively 
dividing at each moment. As each cell 
completes the mitotic cycle, it either differ¬ 
entiates, making a transition into the first 
CFU-GM proliforative compartment, or it 
passes immediately into Gl, starting the 
next round of cell division, or it lapses into 
GO, the resring state. The fraction that dif¬ 
ferentiates is determined by the current 
numbers of resting and cycling stem cells. 

CFU_binh_rate- fractional birth rate of 
CFU-GM cells (average birth rate pet CFU- 
GM cell per unit time). If the more mature 
GM lineage downstream from these GM- 
comrnitted stem cells becomes depleted, the 
division (i.e., birth) rate of the CFU-GM 
cells increases to compensate. Conversely, if 
the downstream cell populations are well 
stocked, then the birth rate among CFU- 
GM cells slows. Tile CFU-GM birth rate 
is also affected by feed-forward from the 
earlier stem cell compartment. 

P_binh_raie = fractional birth rate of 
proliferative cells. These cells are assumed 
to respond only to signals from CFU-GM 
and later cells, up through and including 
the bone marrow reserve, but not to signals 
from the earlier (pre-CFU-GM) stem cells. 

Release_rate__to_blaod = fractional release 
rate of mature ^anulocytes from bone mar¬ 
row reserve to peripheral blood (per reserve 
cell per unit time). This race is controlled by 
a direct feedback loop from the peripheral 
blood compartment; there is no longer range 
control from upstream cell populations. 

The nonlinear feedback control equa¬ 
tions determining these quantities are sum¬ 
marized in the appendix. The specific 
functional forms used are empirical approx¬ 
imations obtained by fitting smooth curves 
to experimental data and to extreme or 
limiting conditions (e.g., maxima] birth 
rates obtained under experimental condi¬ 
tions of maximum stimulation). Sensitivity 
analysis reveals that the exact forms of the 
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equations are less important than their 
overall shapes. Steinbach et al. (22) discuss 
in greater detail the experimental basis 
and the assumptions used to derive the 
feedback control equations. 

Using the Decomposition 
Framework for Risk Extrapolation 

A principle advantage of the causal decom¬ 
position approach to dose-response model¬ 
ing is that it provides an alternative to dose 
metrics as a basis for risk modeling. Thus, 
even though a prediedvely useful dose met¬ 
ric may not exist for benzene, risks can still 
be quantified and extrapolated. The princi¬ 
ple of extrapolation is illustrated by the 
following formula; 

(j»ll*)b uman, inhalation ( P ^ jO mouse, gavage 

u.man, inhalation 

This formula shows that the relation 
ip\\y) estimated from experimental data 
involving mice exposed via oral gavage can 
be composed with the relation (yllx) esti¬ 
mated from data (or from a PBPK model) 

for humans exposed via inhalation. The 
result is an extrapolated (across species and 
route) estimate of the dose—time—response 
relation {p\\x). Each term in the decompo¬ 
sition is subscripted to indicate die type of 
data from which it has been estimated. 
Ideally, all terms would be estimated from 
the same type of data (e.g., sex, age, spe¬ 
cies, strain, route of administration, etc.), 
namely, the data type to which the risk 
model as a whole is co be applied. In prac¬ 
tice, however, if data for humans exposed 
via inhalation (for example) are not avail¬ 
able to quantify some of the components, 
then animal data may be substituted. This 
provides a method of incorporating human 
data where it is available and using only ani¬ 
mal data where necessary. Models quantified 
in this Way can subsequendy be refined when 
more relevant data become available. 

Results 

Results of Monte-Carlo 
Uncertainty Analysis: The Slope 
of the Dose -Response Curve Is 
Unlikely to Be Positive at the Origin 

Figures 2 to 6 show the results of applying 
the maximum-entropy Bayesian Monte- 
Carlo uncertainty analysis technique to the 
data in Table 1. The posterior cumulative 
distribution functions (CDFs) for the tumor 
response probabilities after conditioning on 


the sample sizes and observed tumor inci¬ 
dence rates are shown in Figure 2. These 
provide the input to the rest of die compu¬ 
tations summarized in Figure 1. Figure 3 
shows the CDF for q x . Notice that there is a 
high probability (about 90%) that q x is neg¬ 
ative, suggesting that there is no excess risk 
at low doses. Figure 4 shows this explicitly: 
as administered dose increases from I 
mg/kg/day to 5 mg/kg/day, the PDF for 
risk shifts left. However, at 25 mg/kg/day, 
the lowest dose group for which data are 
available in Table 1, the risk is significandy 
increased compared to die zero-dose risk. 
The interpretation of negative risk values at 
lower doses Is discussed below. Figure 5 
shows the expected value of the dose- 
response function (the solid curve in the 
middle) with 75 and 95% probability bands 
around ic. These bands are to be interpreted 
in terms of die Monce-Carlo-derived poster¬ 
ior PDF or CDF for the response probabil¬ 
ity at each dose level rather than in terms of 
classical hypothesis testing. Thus, they 
replace the 75 and 95% upper and lower 
confidence limits that would be derived 
from the asymptotic X 2 distribution using a 
traditional q x * approach. Figure 6 shows 
analogous results for excess risk, defined as 
total risk minus background (zero-dose) 
risk. The 95% upper probability band for 
excess risk has a positive slope of 0.004 and 
is linear at low doses. This is about half the 
value of q x * derived using the traditional 
(GLOBAL ’82) approach. 

All these results are based on simple 
polynomial regression as a curve-fitting tech¬ 
nique, as described in the methods section. 
More sophisticated nonparametric smooth¬ 
ing (e.g., splines) could be used instead. 

Results of Traditional Multistage 
Modeling Using PBPK Internal 
Doses: The Dose-Response Curve 
Has Zero Slope at the Origin 

The strategy of decomposing the benzene 
dose-response relation as 

{p\\x) = (p\\y)*{y\\x) 

has been implemented for the data in 
Table 1 (25). The internal dose microrela¬ 
tion (yllx) for mice was quantified in 
three steps: a) The PBPK model of Travis 
et al. (17) was used to calculate the steady- 
state amounts of total benzene metabolites 
produced per day from different amounts 
of benzene administered per day. b) A 
nonlinear (Michaells-Menten) regression 
model was used to quantify the relation 
between x, the administered daily dose of 






Figure 2. Beta posterior distributions for the true 
response probability in each doss group. 



Figure 3. Cumulative probability density function (CDF) 
for (Jt. The CDf is obtained via Monte-Carlo simulation, 
by propagating numbers sampled from the beta distrib¬ 
utions on the left of Figure f through the rest of the 
diagram. The CDF for each derived quantity in Figure 1 
is obtained via an algebraic combination of the (sam¬ 
pled) values of the quantities pointing into it. 
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Figure 4. Sensitivity analysis ot the probability density 
function of risk to different administered concentrations. 
At doses x=0.2 and x=1 mg/kg/day, the probability 
mass is concentrated near zero (and is less than 0.1). 
For x=5 mg/kg/day. the probability mass is shifted 
left: the probability of a positive tumor response is less 
than at lower doses. At x-25 mg/kg/day, however, it is 
very likely that the risk is positive. 


benzene, and y, the total amount of ben¬ 
zene metabolites formed per day. c) The 
fractions of total benzene metabolites 
formed along different specific metabolic 
paths (e.g., leading to hydroquinone conju¬ 
gates, muconic acid, etc.) were quantified 
based on empirical data. 

The mouse PBPK model relation 
between administered dose and internal dose 
(defined as total benzene metabolites formed 
per day mg/kg/day in mice exposed via oral 
gavage to X mg/kg/day of benzene) was well 
approximated by the following nonlinear 
(Michaelis-Menten) regression model: 

O' 1 *) mount, gavage = 76.4x/(80 .75 + x) 

(x in mg/kg/day via gavage) 

The corresponding relation for humans 
exposed to benzene inhalation derived 
from a benzene PBPK model for humans 
(17) was 

O' 1 x) human, inhalation = 9.525*/(I 32.24- x) 
(xin ppm via inhalation). 

Orher choices of route and dose regimen 
lead to similar formulas. Each such formula 
gives a useful numerical approximation of 
the results of the full dynamic PBPK 
model for a specific choice of che exposure 
summary variable, x, the metabolism vari¬ 
able, y, and for a range of administered 
dose scenarios. 

Next, the microrelarion (filly) was 
quantified using a traditional multistage 
dose-response model, talcing both total 
benzene metabolites and the amounts of 


Figure 5. Probability bands for the doss-response 
relation. For doses below about 1C mg/kg/day, it is 
unlikely that there is a positive excess lumor risk due 
to benzene exposure. 


metabolites formed along specific meta¬ 
bolic pathways (corresponding to specific 
components of y) as the dose variable, and 
using the tumor data in Table 1 for 
response data. The internal dose versus 
observed response data points for total 
benzene metabolites are shown in Table 3. 

The projected internal dose levels (the 
y* values) in Table 3 are estimated from 
the corresponding administered dose levels 
(x = 0, 25, 50, and 75 mg/kg/day) using 
the Michaelis-Menten statistical interpola¬ 
tion model y*(x) = 76.4x/(80.75 +x). The 
traditional multistage methodology (imple¬ 
mented in GLOBAL ‘82) applied to the 
internal dose estimates in Table 3 produced 
the following results: 

The MLE Cancer probability from y 
mg/kg/day of benzene metabolites (for 
male B6C3F1 mice) is 

(p\\y)~ 1 -exp(-0.001)013j 3 ) 

(MLE internal dose—response 
model for mice). 

The estimated 95% UCL for die low-dose 
slope, calculated via the q* mechodology, is: 

4, *=0.00554. 

(estimated 95% UCL). 


Table 3. Estimated internal dose-response functions 
(or mice. 


Administered dose x 
mg/kg/day 

0 

75 

50 

75 

Estimated dose. 

/ Img/kg/day) 

0 

IB 

29 

42 

% of animals with 

0 

8 

40 

74 


squamous cell carcinomas 
(50 animals/dose group) 


Figure 6. Probability bands for excess risk. The 
expected excess risk (solid line) is negative at doses 
below 15 mg/kg/day. However, there is a 5% (but not 
a 25%) chance that it is positive even at very low 
doses, as predicted by the qT methodplegy. 


For humans, the dose-response function is 
estimated as follows: 

(pi lx)h uman, inhalation 'P ^/human, inhalation 
(7 ^ *0 human, inhalation 
(fi 11 "/)human, inhalation 

*9.525x/(132.2+x). 

In the absence of specific data from which 
to estimate (^>l!.>) hllmi „, motion- the inter¬ 
nal dose—response function estimated for 
mice (the most sensitive known species, 
sex, and strain) was used instead: 

estimated (jpll Jf)hum ln , inhalation 
= (p ll.kLoust. gavage* 9.525x/(l 32.2 + *) 

= 1 —exp{—0.000019[9-525-*/(I32.2 + .v)] 3 } 
(estimated human dose-response function). 

Note that the MLE model for 
(/■lljkOmomt, gavage has a slope of zero as y 
approaches zero (and hence as x approaches 
zero). The estimated total dose-response 
relation (p\! ar/human, inhalation is nonlinear 
with a slope of zero at x = 0. A maximum 
entropy Bayesian Monte-Carlo uncertainty 
analysis of the mouse internal dose- 
response function (p Ily) mouse , Eavige con¬ 
firmed that there is no increase in expected 
risk at doses below about 14 mg/kg/day. 

Applications c»f the Human Risk 
Model to Different Dosing Scenarios 

Tables 4 and 5 show the results of applying 
the foregoing risk model to two human 
exposure scenarios of practical interest: a 
daily dosing scenario (Table 4) and a contin¬ 
uous exposure scenario (Table 5). In Table 
5, the risk to humans from routine expo¬ 
sure is shown as a function of an unknown 
constant, k. This is defined as the ratio of 
the steady-state, biologically effective dose 
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Table Estimated human dose-response model for 8 hr/day lifetime inhalation exposures to benzene. A 



B 


Table 5. Individual excess cancer risks from lifetime inhalation exposure to 1 ppm of benzene, for different values 5.95e+oa 
of Xr(= mouseihuman internal dose ratio) 


Doss measure 

Risk model 


Value of k 


Z.98e+08- 

f\ 

0.01 

0.1 

1 

L 

Administered dose 

1.45x10 4*+1.31 *10-4*2 

D 0016 

0.0016 

0.001 B 



Total metabolites 

1.86 xKTV 

1.1 E-13 

1.1E-10 

1.1 E-7 

| 

i\\ — 

Total HQC dose 

1.51x1(riV 

0.8 E-13 

0.8 £-10 

O.B E-7 


l\ v - - 

Total MA dose 

1 51 x1C- J y 3 

8.1 E-13 

8.1 E-IQ 

8.1 E-7 

0.00 - 1 

- -- --— r W T — ■V T» | 

Total metabolites, two-stage IMS 

6.21 x 10-4 yl 

2.0 E-9 

2.0 E-7 

2.0 E 5 

0.00 240.00 400.00 720.00 98D.Q0 









Abbreviations: HQC, hydrnqurnone conjugates; MA, muconic acid; LMS, linearized multistage risk model. 


Table 6. Maximum likelihood linearized multistage risk models fit to male mouse squamous cell cancer incidence 
data for different surrogate internal dose measures andiumor end points. 


y= Dose surrogate 

Maximum likelihood LMS risk model 

Description 

End point—all squamous cel! carcinomas 

Administered dose 

1.45x 10 -3 *+ 1.31 x 1(H 

Linear-quadratic 

Total metabolites 

1.B6x KT 5 y 3 

Cubic 

Total HQC dose 

1.51 xKrV 

Cubic 

Total MA dose 

1.51 x 10" 2 

Cubic 

AUC HQC dose 

4.14 y z + 2.44 y 3 

Quadratic-cubic 

AUC MA dose 

0.448 y 2 +0.004 y 3 

Quadratic-cubic 

End point—squamous cell carcinomas 
of Zymbal and preputial glands 3 

Zymbal gland SCC 

6.13xirV 

Cubic 

Total metabolites 

Preputial gfand SCC, total metabolites 

1.51 x 1<Hy 3 + 8.89 x hr 6 y 3 

Quadratic-cubic 


Abbreviations. IMS, linear multistage risk models; SCC, squamous cell carcinoma. "Specific tumor type end points 
should only be analyzed by time-to-tumor models that adjust for competing risks due to other tumors. These quanta! 
models are for purposes of illustration only. From Cox and Ricci 1 20). 


of the carcinogenic metabolite in humans 
to its value in mice, assuming low, constant 
administered concentrations of benzene. 
Tabic 6 shows that the main finding of a 
zero slope for the interna] dose—response 
function at the origin (i.e„ low-dose non¬ 
linearity) holds for many choices of inter¬ 
nal doses and end point. Details of the 
computations summarized in Table 6, are 
presented by Cox and Ricci (20). 

Results of Computer Simulation 
Modeling: the Hcmatotosicity 
Dose-Time-Response Relation Has a 
Negative Slope at Low Concentrations 

Ideally, the results of risk calculations 
based on the decomposition 

(p\\ x) = (p\\y) * (jrllx) 


would be cross-checked and improved using 
the results of alternative decompositions 
such as 

(p liar) = (pll IV, p) * (iVjplIx). 

Although some data are available to quantify 
the cytotoxic effects on Wand the cytotoxic 
and genotoxic effects on p of inhalation 
exposure to benzene (7), a detailed calcula¬ 
tion of benzene dose-response relations by 
this decomposition analogous to the one 
based on PBPK modeling of internal doses, 
has not been done. Nonetheless, it is now 
possible to provide a key component of 
such a model by using a simulation model 
of benzene-induced hematotoxicity. 
Figures 7 to 3 show results from the bio¬ 
logically based computer simulation model 


. 10 ppm berrcene 

-- - 25 ppm benzene 

- 50 ppm benzene 

-1DQ ppm benzene 

-- 200 ppm bepzene 

Figure 7. Computer simulations of early and lata 
CFU-GM responses to different benzene concentra¬ 
tions, Curves show the simulated time courses of 
CFU-GM cells—Ml earliest CFU-GM 1 : (fl) most 
mature CFU-GM 10—in response to constant expo¬ 
sures to (approximately) ID, 25, 5Q, 1D0, and 200 ppm 
of benzene, respectively. 



. 24 hr between repealed doses 

--72 hr between repeated doses 

- 144 hr between repeated doses 

Figure 8. Different amounts of spacing between 
repeated 8-hr inhalation doses of benzene cause very 
different hematopoietic responses for the same admin¬ 
istered dose (AUC| per unit time. Curves correspond to 
24, 72, and 144 hr between repeated doses. 

used to quantify the relation (AMI*). There 
are two major conclusions from this 
model, a) The function N(x) = (Ad|jv), 
expressing the steady-state size of early 
myeloid proliferative cell population (e.g., 
HPCs or CFU-GM) as a function of the 
administered benzene concentration 
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Figure 9. Steady-state relations between benzene exposure and hematopoietic ceil population sizes predicted by 
a dynamic simulation model. Benzene doss axis may need to be rescaled for humans. 


assuming a constant exposure scenario, is a 
nonmonotonic function of benzene con¬ 
centration (Figure 9). At sufficiently low 
concentrations, benzene exposure is pre¬ 
dicted to increase cell population sizes (for 
early HI*Cs and immature CFU—GM cells, 
but not for very mature CFU—GM cells). 
In Figures 7 to 9, CFU-GM1 is the earli¬ 
est of the 10 age-structured subcompart¬ 
ments of the CFU-GM population, while 
CFU-GM10 is the Iasc of these subcom¬ 
partments. Figure 9 plots the steady-state 
levels in Figure 7 against administered con¬ 
centration. It shows that che relation 
between ppm of benzene and predicted 
steady-state values of N (where Nis one of 
the early CFU-GM or HPC populations) 
is positive at sufficiently low concentrations 
but negative at higher concentrations. If it 
is conjectured that AML risk is increased 
only by exposure scenarios that signifi¬ 
cantly depress CFU-GM and/or HPC 
populations (e.g-, because recruitment into 
active cycling or clonal expansion of pre¬ 
leukemic stem cells becomes more likely), 
then the predicted shape of the {/Vlljr) 
relation would be consistent with the 
hypothesis that benzene at sufficiently low 
concentrations does not increase (and 
might even decrease) AML risk, b) The 
same AUC of administered dose may have 
profoundly different effects on hematotoxi- 
city depending on how it is administered 
over time. For the same AUC per unit time 
of administered benzene, higher concentra¬ 
tions may lead to disproportionately large 
hematotoxic effects (Figure 8). The 


nominal administered concentration in 
Figure 8 is 100 ppm. Different concentra¬ 
tions wirh the same time pattern of admin¬ 
istration produce similar results for a broad 

range of administered concentrations. Since 
the benzene hematotoxic potency factors in 
the model were originally developed for 
mice rather than for humans, however, 
human data are needed to calibrate the 
model more exactly for benzene. 

Although the simulation model pro¬ 
duces a wealth of other results, it should be 
developed and validated further before 
being used for risk assessment. Nonethe¬ 
less, the model-based conclusions appear to 
be consistent with the results arrived at by 
other methods. 

Discussion 

Bayesian Monte-Carlo Uncertainty 
Analysis Results 

The 95% upper probability bands on low- 
dose slopes for excess risk derived by 
Bayesian Monte-Carlo uncertainty analysis 
(Figure 6) are of the same order of magni¬ 
tude as the (j\ * values derived by traditional 
methods. However, the Monte-Carlo 
uncertainty analysis in Figures 3 to 6 also 
reveals a great deal of additional informa¬ 
tion that many decision makers—includ¬ 
ing all who follow the expected utility 
paradigm of rational decision making 
(23 )—would consider relevant. Perhaps 
most striking is the fact that the expected 
value of the excess risk is negative for 
administered doses smaller than about 


15 mg/kg/day (Figure 6). Taken at face 
value, this finding suggests that the health 
protection benefits from reducing doses 
below this level are very uncertain, and may 
not be significantly different from zero. 
Notice that while the total risk of tumor 
defined as a probability or hazard rate, must 
be nonnegative, the risk due to a specific 
cause such as benzene can perfectly well be 
negative. For example, risk might be 
reduced if benzene exposures in the 1 to 15 
mg/g/day range were to suppress the 
expression of carcinogenic damage (e.g., by 
killing or preventing the division of affected 
cells) or were to reduce the sizes of the pop¬ 
ulations at risk of carcinogenic damage. 
Since benzene was once used as a human 
chemotherapeutic agent, such effects on cell 
kinetics would not be completely unex¬ 
pected. Without speculating about such 
mechanisms, the interpretation of the nega¬ 
tive risk values in Figures 4 through 6 is 
that the expected lifetime tumor probability 
at low doses is less than the expected life¬ 
time rumor probability at zero dose. 

Traditional Risk Model 
Based on Internal Doses 

The benzene risk model based on PBPK 
modeling and use of internal doses, namely, 

(PI Inhuman, inhalation 

1 -expf-0.000019[9.525x/(132.2+*)] 3 l, 

differs significantly at low dose from the 
traditional {administered dose) benzene 
risk model discussed in the introduction. 
For example, at x= 12 mg/kg/day, mice 
form approximately y= 76.4* 12/(80.75 + 
12)= 10 mg/kg/day of total benzene 
metabolites, giving a projected lifetime 
cancer risk of , p *(10) - 0.019 based on the 
internal dose model. By contrast, the 
administered dose risk model p*(x) = 

1 — expf—(0.00145* + 0.00013x- 2 ) gives a 
predicted lifetime cancer risk level of 
p*(12) =0.035 based on the administered 
dose model. Thus, fitting the dose-response 
function to internal doses instead of admin¬ 
istered doses changes the maximum likeli¬ 
hood linearized multistage model from 
linear-quadratic to cubic, reducing estimated 
risk by about a factor of 2 at an administered 
dose level of x= 12 mg/kg/day. [This oral, 
gavage dose level corresponds roughly to a 
6-hr inhalation dose of about 10 ppm based, 
on total metabolites formed (L9).] Similar 
calculations for x = 1 mg/kg/day give 
p*(0. 93) = 0.000015 for the Michaells- 
Mentcn model for yi*(l) = 0.000019 for 
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the PBPK model] compared to />*(1) = 
0.0016 for the administered dose model, 
i.e., at an administered dose level of 1 
mg/kg, risk estimated from the internal 
dose model is about 1% of the risk esti¬ 
mated from the administered dose model. 
At lower doses (corresponding to inhalation 
doses below 1 ppm) the cubic internal dose 
model gives lifetime cancer risk predictions 
less than 1% as large as those from the 
linear-quadratic administered dose model. 

There are several reasons to believe that 
the cubic interna! dose—response risk model 
may be more realistic than the adminis¬ 
tered dose—response, model. If a Weibull 
model of the form J>{y) = I — exp[—(m by c )] 
is St to the four internal dose response 
points in Table 4, the result is the MLE 
Weibull risk model 

(/>lly) = 1 -exp(- 0.000015/ 3065 ). 

Thus, the conclusion that the relation 
between interna! dose and probability of 
response is approximately cubic (rather 
than linear or quadratic) does not appear 
to be only an artifact of using the multi¬ 
stage family of risk models. Also, the inter¬ 
nal dose risk model provides a slightly 
better fit than the administered dose risk 
model to the available animal dosc- 
response data. The sample values of (p II x) 
at x = 25, 50, and 100 mg/kg/day are 
(/>llx) = 0.08, 0.40, and 0.74. The corre¬ 
sponding values predicted by the adminis¬ 
tered dose linear-quadratic model are 0.11, 
0.33, and 0.77, while the values predicted 
by the cubic model with the Michaelis- 
Menten internal dose model are 0.105, 
0.37, and 0.755. In each Case, the best¬ 
fitting (cubic) internal dose—response 
model predicts values of (pllx) slightly 
closer to the observed ones than the corre¬ 
sponding predictions from the best fitting 
(linear-quadratic) administered dosc- 
response model. Quantitatively, the inter¬ 
nal dose model provides a better fit to the 
data, with the maximized value of the log- 
likelihood function being -76.53 for the 
internal dose model compared to —77.15 
for the administered dose model. Finally, 
the internal dose-response model appears 
to be much more consistent with epi¬ 
demiological data on rumors in heavily 
benzene-exposed workers (Turkish shoe 
workers and Chinese workers) than the 
administered dose—response model. 


For an administered dose of x = 1 
mg/kg/day, the internal dose model 
predicts an MLE risk of oniy 0.0000155, 
which is about 1% as large as the MLE risk 
predicted by the administered dose model. 
On the other hand, the respective upper 
UCLs for risk ai x - 1 mg/kg/day are 0.005 
for the internal dose model compared to 
0.008 for the administered dose model. 

Thus, even though the MLE estimates dif¬ 
fer by a factor of 100, the 95% UCLs 
based on traditional q\ uncertainty analy¬ 
sis are very similar. This insensitivity of 
95% UCLs to significant variations in risk 
models is typical; it has been both criti¬ 
cized as a defect and praised as a virtue of 
the ij\ * approach. 

Biologically Based 
Risk Assessment Modeling 

Biologically based computer simulation 
models are available for quantifying some 
of the key micro relations thought to be 
involved ill benzene leukemia causation, 
including (yll x) and (AMI x). Where such 
models exist, they can be used to quantify 
the output histories, like [y\ and {IV], 
caused by any benzene administered dose 
history, {xj. Instead of using dose metrics 
for risk extrapolation, therefore, it is 
possible to construct explicit computer 
simulation models of the pharmacokinetic 
and pharmacodynamic processes (ideally 
including hematotoxicity and genotoxicity) 
mediating between benzene exposure and 
AML induction. 

Such biologically based risk assessment 
(BBRA) models have the advantage of mak¬ 
ing explicit, testable predictions of observ¬ 
able quantities (such as quantities of 
benzene metabolites formed, cell population 
sizes, and so forth) as well as tumor risk, 
which is much harder to observe directly. 
They also bypass the need to rely upon dose 
metrics, which may not exist in any useful 
form for benzene. On the other hand, 
BBRA models are often too complex and 
detailed and require too many uncertain 
inputs to be used for regulatory risk assess¬ 
ment and public health risk management. 

A principle expoiced in this article to 
avoid complexity is that if administered 
dose histories are restricted to simple 
forms, then the full dynamic simulation 
approach becomes unnecessary. Specifi¬ 
cally, under a wide range of conditions, a 
constant administered dose eventually 
produces a constant steady-state output for 


each component of ly) and {/V{. In this 
case, letting X denote the constant applied 
concentration (in ppm, for inhalation 
exposures) and letting y(x) and N(x) 
denote the resulting steady-state levels of a 
specified benzene metabolite and cell pop¬ 
ulation, respectively, it is possible to quan¬ 
tify the relations (y II x) and (AMI x) using 
nonlinear regression models such as 

y(x) = Vx/{K+ x ) 
(Michaelis-Menten nonlinear 
regression model) 

where V and K are model parameters to be 
estimated from observed data (such as 
[x,y\ pairs, in this example). Alternatively, 
if a PBPK model is available, as is the case 
for benzene, it can be used to quantify 
many ( x , y ) pairs, and then a relatively 
simple nonlinear regression model can be 
used to describe rhe results, bypassing the 
need to use the PBPK model itself when 
only simple, constant dose scenarios are 
being considered. 

Conclusion 

The conclusions of greatest potential prac¬ 
tical interest from the preceding analyses 
are as follows, a) There is no evidence of a 
positive relation between benzene exposure 

and tumor probability (or causal antece¬ 
dents such as depressed myeloid stem cell 
and HPC populations) at benzene concen¬ 
trations below 1 ppm. b) Different analytic 
approaches — including Bayesian Monte- 
Carlo uncertainty analysis and PBPK-based 
internal dose modeling—suggest that the 
(/’ll x) curve relating benzene concentra¬ 
tion to AML risk at sufficiently low, con¬ 
stant concentrations of benzene approaches 
a zero or negative slope as concentration 
falls below abour 10 ppm. c) Empirical 
data and biologically based risk models 
agree thac, for the same total administered 
dose (AUC), higher concentrations of ben¬ 
zene may cause disproportionately large 
hematotoxic responses. 

Quantitatively, none of the statistical or 
biological evidence investigated suggests 
that reducing benzene concentrations 
below 1 ppm would have any detectable 
health benefits. The hypothesis of a zero or 
nonpositive slope for excess risks due to 
inhalation of low concentrations of ben¬ 
zene may have sufficiently interesting 
public health policy implications to merit 
further evaluation. 
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Appendix: Benzene 
Hematotoxicity Model Equations 
and Parameter Values 

HPC Cells 

d (cycIing_HPCs)/«/r = recruited_from_rest!ng + (cell_divisions - 
HPC_deaths) - return_to_resting - differentiated 

d fresting_HPCs)/V/ = return_to__resring - recruited_from_resting 
recruited_from_resting = HPC_reeruitment_rate * 
resting MFCs 

cclLdivisions = recruired_from_re$ting {since each recruited 
Cell divides at mitosis} 

111’C_ deaths = cycling_HPCs * cytotoxic_potency„for_HPCs 
* metabolire__in_ce!ls 

metabolite_in_cells = output from PK mode! 

rerurn_to_resting = f 1 - HPC_<Ufferentiati<m_fractioti) * 
cycling_HPCs / HPC_cycle_time 

differentiated = HPC_differentiation_fraction * cyding_HPCs/ 
HPC_cyde_time 

CFU-GM Cells 

r/CFU(l, t)/dt = differentiated-[CFU(1, r)/CFU_transit_time] - 
(CFU_death_rate - CFU_birch_rate) * CFU(1, t) 

dCk\J(i,t)!dt = [CFU(r-1, r)/CFU_transit_time] - 
[CFU(r, t )/(iFU_ transit_time! - (CFU_deach_rate - CFU_birth_ 
rate) * CFU(4 t ), for i = 2, 3,10 

GM Proliferation 

dP(\,t) = inflow_from_CFU_GM(t) - [/ 3 (l ) r)/ , / , _transit_time] - 
(/Cdeuth. race /“ birth rate) * P(l,t ) 

dP(i,t) = [P(i — l,f)//Ltransit_timeJ — [P(i, f)//Ctransit_rime] - 
(/“ dcatli .rate - / J _birth rate) * P(i, t), for j = 2, 3, .... 10 

inflow_from_CFU-GM(r) = G(10, t)/CFU_transit_time 

/ J death_ rate = cytotoxidty_potency_for_proliferacing_celis * 
metaboli ce_i n_cells 

Nonproliferative Populations 

dMIdt- prccursors_in - tmmred_ceIls_out 

precursors_in = matured_ceIIs_out(f) = M(r)/macutation_tran- 
sit_time 

tiRldt = matured_cells_out - release_rate_to_blood * R 
dB/dt- release_rate_to_blood * R— decay_rate*5 


HPC_recruitment rate = 

k\ * EXP(-vl * resting_HPCs — A * cyding_HPCs) + k2 * 
EXP{—k2 * 7.) h 

*1 = 0.36 {per hr} kl - 0.04 {per hr} 

kb = 0.008 {per hr} 

v\ - 1.6955E-6 vl = 3.1918E-11 

Z = {signal_Z + CFU_GM) {=3.208E9 + 3.469E10 + 
12.8E10 in sceady state} 

signal_Z - proliferating_pooI + maturation_pool + mar- 
row_reserve = P+ M+ R 

HPC diffcrelltiation fraction - 

1 — MAX(0, 0.5 + 1.76E-9 * (sceady _s tat e_HPCs — rest- 
ing_HPCs - cyding_HPCs)) 

steady_state_HPCs = 3.125E8 {HPC cdls in sceady state} 

CFU_birth_rate = 0.1/(3 + min(0.2, kb * testing_HPCs + k4 * 
cycling_HPCs) +- min(0.2, kb * CFU_GM) + min(0.2, kG * 
signa1_Z)) 

kb = 3.6418£-9 kb = 3.5476E-10 

kG = 6.9953E-12 

(22, calculated from steady-state conditions) 

CFU_GM = CFU_GM_ 1 + CFU_GM_2 + CFU_GM_3 + 
CFU_GM_4 + CEU_GM_5 + CFU_GM_6 + 

CFU_GM_7 + CFU_GM_8 + CFU_GM 9 + 

CFU_GM_I0 

CFU_transit_time =10 {hr for each of 10 compartments = 
100 hr} 

CFU_GM_cytotoxicity_death_rate = metabolite_in_ccIIs * 
ctroroxic._poteucy_for_CFU.GiVl 

cytotoxic_potency_for_CFU_GM = 0.13 {per proliferating 
cell/hi} 

/Murth rate - 0.07/(l + min(0.6, k7*Z)) (19, eq. [16]) 
kl = 3.7289666E-12; Z= (signal_Z+ CFU_GM) 
/{_transir_time =5.6 [hr for each of 10 compartments = 56 hr} 
cytotoxic_potency_for_proliferat:ing cells = 0.13 {per prolifer¬ 
ating cell/hr} 

P = prol iferating_pooI = PI + P2 + P3 + P4 + P5 + P6 + P7 
+ P8 + P9 + P10 

/*_death_tate = cytotoxkcy_death_rate_for_pro!iferating_cells = 
metabolite in ceils * cytotoxic_potency_for_prohferacing 
cytotoxic_potency_for_proliferating =0.13 {per prolifer¬ 
ating cell/hr} 

release_rate_to_hfood = kS - kb> * EXPI-i’3 * GRF) = 2 — 
1.9875 * EXP(-6.309E-4 * GRF) 
blood_decay_rare = 0.1 (22) 

GRF_decay_rate = 0.1 {per hr (22)) 
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